library(data.table)
library(estimatr)
library(fixest)
library(marginaleffects)
library(ggplot2)
library(ggforce)
library(lubridate)
library(interflex)
library(tinytable)
library(texreg)
library(xtable)

library(here)

here::i_am("scripts/analysis.R")

# Functions ---------------------------------------------------------------

get_cis <- function(model, params, levels = c(.95, .90), lh = FALSE) {
  
  if (lh) {
    cis <- data.frame(coef = names(params),
                      est = summary(model)[params, "Estimate"],
                      se  = summary(model)[params, "Std. Error"])
  } else {
    cis <- data.frame(coef = params,
                    est = coeftable(model)[params, "Estimate"])
  }
  
  
  for (ci in levels) {
    cis <- cbind(cis, confint(model, 
                              parm = params,
                              level = ci)  )
  }
  
  cis$var_name <- rownames(cis)
  return(cis)
}

# Figure 1 -------
#https://fred.stlouisfed.org/series/WCBLSA

data <- fread(here('data', 'WCBLSA.csv'))

data$amount <- data$WCBLSA / 1000

ggplot(data[year(DATE) >= 2005 & year(DATE) < 2022],
       aes(x = DATE, y = amount)) +
  geom_ribbon(aes(ymax = amount, ymin=0), 
              color = 'black', fill = 'blue', alpha = .6) +
  scale_y_continuous(expand = expansion(mult = c(0.002, .03))) +
  scale_x_date(date_breaks = '2 years',
               date_labels = "%Y") +
  theme_bw(base_size = 14) +
  labs(x = element_blank(),
       y = '$USD (billions)') +
  facet_zoom(xy = year(DATE) > 2014 & year(DATE) < 2020,
             horizontal = FALSE,
             zoom.size = .75) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        zoom.y = element_rect(fill = NA, color = NA))

ggsave(here('output', 'figures', 'figure-1.pdf'),
       height = 4, width = 7, units = 'in')
rm(data)

# Load data ---------------------------------------------------------------

survey1 <- fread(here("data", "survey-1-cleaned.csv"))

survey1$fed_treatment_policy <- factor(survey1$fed_treatment_policy,
                                       levels = c("fed_treat_lend",
                                                  "fed_treat_dem",
                                                  "fed_treat_mh",
                                                  "fed_treat_risk"))

survey2 <- fread(here("data", "survey-2-cleaned.csv"))


# Table 1 -----------------------------------------------------------------

treats <- table(survey1$fed_treatment_china[survey1$attentive == 1], 
                survey1$fed_treatment_policy[survey1$attentive == 1])

treats <- cbind(treats, apply(treats, 1, sum))
treats <- rbind(treats, apply(treats, 2, sum))

printTable <- xtable(treats)

print(printTable, file = here("output", "tables", "table-1.tex"))

# Figure 2 ----------------------------------------------------------------

## Analysis ---------------------------------------------------------------

models_svy1 <- list()

### Policy Support, Binary, Unweighted -----
svy1_policy_bin <- lh_robust(
  outcome_policy_support ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi,
  data = survey1[survey1$attentive == 1, ],
  linear_hypothesis = c("fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh",
                        "fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk"),
)

models_svy1 <- append(models_svy1, list("Support, Binary, Unweighted" = svy1_policy_bin))

### Policy Support, Binary, Weighted -----
svy1_policy_bin_w <- lh_robust(
  outcome_policy_support ~ fed_treatment_policy + fed_treatment_china 
  + democrat + republican + college + female + stock + age + hhi,
  weights = w_ebal_attentive,
  data = survey1[survey1$attentive == 1, ],
  linear_hypothesis = c("fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh",
                        "fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk"),
)

models_svy1 <- append(models_svy1, list("Support, Binary, Weighted" = svy1_policy_bin_w))

### Policy Support, Integer, Unweighted -----
svy1_policy_int <- lh_robust(
  outcome_policy_int ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi,
  data = survey1[survey1$attentive == 1, ],
  linear_hypothesis = c("fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh",
                        "fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk"),
)

models_svy1 <- append(models_svy1, list("Support, Integer, Unweighted" = svy1_policy_int))

### Policy Support, Integer, Weighted -----
svy1_policy_int_w <- lh_robust(
  outcome_policy_int ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi,
  data = survey1[survey1$attentive == 1, ],
  weights = w_ebal_attentive,
  linear_hypothesis = c("fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh",
                        "fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk"),
)

models_svy1 <- append(models_svy1, list("Support, Integer, Weighted" = svy1_policy_int_w))

### Trust, Integer, Unweighted -----
svy1_trust <- lh_robust(
  outcome_trust ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi,
  data = survey1[survey1$attentive == 1, ],
  linear_hypothesis = c("fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh",
                        "fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk"),
)

models_svy1 <- append(models_svy1, list("Trust, Integer, Unweighted" = svy1_trust))

### Trust, Integer, Weighted -----
svy1_trust_w <- lh_robust(
  outcome_trust ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi,
  weights = w_ebal_attentive,
  data = survey1[survey1$attentive == 1, ],
  linear_hypothesis = c("fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh",
                        "fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk"),
)

models_svy1 <- append(models_svy1, list("Trust, Integer, Weighted" = svy1_trust_w))

## Prolific Sample -------

models_svy2 <- list()

svy2_policy <- lm_robust(policy_support_bin ~ treat 
                         + party 
                         + female
                         + age
                         + stock 
                         + income_int
                         + college_bin, 
                         survey2[survey2$attentive == 1,])

models_svy2 <- append(models_svy2, list("Support, Binary, Unweighted, Prolific" = svy2_policy))

svy2_policy_w <- lm_robust(policy_support_bin ~ treat 
                         + party 
                         + female
                         + age
                         + stock 
                         + income_int
                         + college_bin, 
                         weights = w_ebal_attentive,
                         survey2[survey2$attentive == 1,])

models_svy2 <- append(models_svy2, list("Support, Binary, Weighted, Prolific" = svy2_policy_w))

svy2_policy_int <- lm_robust(policy_support_int ~ treat 
                             + party 
                             + female
                             + age
                             + stock 
                             + income_int
                             + college_bin, 
                             survey2[survey2$attentive == 1,])

models_svy2 <- append(models_svy2, list("Support, Integer, Unweighted, Prolific" = svy2_policy_int))

svy2_policy_int_w <- lm_robust(policy_support_int ~ treat 
                               + party 
                               + female
                               + age
                               + stock 
                               + income_int
                               + college_bin, 
                               weights = w_ebal_attentive,
                               survey2[survey2$attentive == 1,])

models_svy2 <- append(models_svy2, list("Support, Integer, Weighted, Prolific" = svy2_policy_int_w))

svy2_trust_int <- lm_robust(trust ~ treat 
                             + party 
                             + female
                             + age
                             + stock 
                             + income_int
                             + college_bin, 
                             survey2[survey2$attentive == 1,])

models_svy2 <- append(models_svy2, list("Trust, Integer, Unweighted, Prolific" = svy2_trust_int))

svy2_trust_int_w <- lm_robust(trust ~ treat 
                            + party 
                            + female
                            + age
                            + stock 
                            + income_int
                            + college_bin, 
                            weights = w_ebal_attentive,
                            survey2[survey2$attentive == 1,])

models_svy2 <- append(models_svy2, list("Trust, Integer, Weighted, Prolific" = svy2_trust_int_w))


## Plot results -----------------------------------------------------------

params_svy1 <- c('Accountability' = 'fed_treatment_policyfed_treat_dem',
                 'Moral Hazard' = 'fed_treatment_policyfed_treat_mh',
                 'Risk' = 'fed_treatment_policyfed_treat_risk',
                 'China' = 'fed_treatment_china')

plotData <- data.table()
results <- data.table()

for (m in 1:length(models_svy1)) {
  results <- rbind(results,
                   cbind("model" = names(models_svy1)[m],  
                         get_cis(model = models_svy1[[m]]$lm_robust,
                                 params = params_svy1,
                                 levels = c(.95, .90))
                   ) # end cbind
  ) # end rbind
}

for (m in 1:length(models_svy2)) {
  results <- rbind(results,
                   cbind("model" = names(models_svy2)[m],  
                         get_cis(model = models_svy2[[m]],
                                 params = c('China\n(Prolific)' = "treat"),
                                 levels = c(.95, .90))
                   ) # end cbind
  ) # end rbind
}

results$outcome <- sapply(strsplit(results$model, ", "),
                          "[[", 1)
results$outcome_measure <- sapply(strsplit(results$model, ", "),
                                  "[[", 2)
results$weight <- sapply(strsplit(results$model, ", "),
                         "[[", 3)

results$attentive <- ifelse(grepl("inattantive", results$model,fixed = TRUE),
                            "Yes", "No")

results$sample <- ifelse(grepl("Prolific", results$model,fixed = TRUE),
                            "Prolific", "Lucid")

results$var_name <- factor(results$var_name, levels = c('China\n(Prolific)',
                                                        'China',
                                                        'Risk',
                                                        'Moral Hazard',
                                                        'Accountability'))

results$outcome_label <- paste0(results$outcome, " (", results$outcome_measure, ")")

ggplot(results, aes(y = var_name, x = est, 
                    color = weight,
                    shape = weight)) +
  geom_vline(xintercept = 0, color = 'red', linetype = 'dashed', size = .5) +
  geom_point(size = 3.5,
             position=position_dodge(width=-1)) +
  geom_errorbar(aes(xmin = `2.5 %`, xmax = `97.5 %`),
                #width = c(rep(.2, 6), rep(.4, 4)), 
                width = 0,
                size = .66,
                position=position_dodge(width=-1)
  ) +
  geom_errorbar(aes(xmin = `5 %`, xmax = `95 %`),
                width = 0, 
                size = 1.45,
                position=position_dodge(width=-1)
  ) +
  scale_color_manual(name = "Weighted?",
                     values = c('black', 'gray')) +
  scale_shape_manual(name = "Weighted?",
                     values = c('circle', 'triangle')) +
  theme_bw(base_size = 12) +
  labs(y = element_blank(),
       x = element_blank(),
       color = 'Sample') +
  theme(legend.position = c(0.75, 0.2),
        legend.title = element_blank()) +
  facet_wrap(~outcome_label, 
             ncol = 2,
             scales = "free_x")

ggsave(here("output", "figures", "figure-2.pdf"),
       width = 6.5, height = 5, units = "in")


# Table 2 -----------------------------------------------------------------

results <- data.table()

for (m in 1:length(models_svy1)) {
  results <- rbind(results,
                   cbind("model" = names(models_svy1)[m],  
                         "coef"  = c("Moral Hazard", "Risk"),
                         tidy(models_svy1[[m]]$lh)[ , c("estimate", "std.error", "p.value")]
                   ) # end cbind
  ) # end rbind
}

results$outcome <- sapply(strsplit(results$model, ", "),
                          "[[", 1)
results$outcome_measure <- sapply(strsplit(results$model, ", "),
                                  "[[", 2)
results$weight <- sapply(strsplit(results$model, ", "),
                         "[[", 3)


results_w <- results[results$weight == "Weighted", -"model"]
results_uw <- results[results$weight == "Unweighted", -"model"]

table_2 <- t(results_uw[ , c("estimate", "std.error", "p.value")])
table_2 <- round(table_2, 3)

table_2[2, ] <- paste0("(", table_2[2,], ")")

table_2 <- cbind(data.frame("$\\beta_X:$" = c("Estimate", 
                                              "\\textit{SE}", 
                                              "\\textit{p-value}"), 
                            check.names = FALSE),
                 table_2)

colnames(table_2)[2:7] <- paste0("$\\beta_{\\text{", results_w$coef, "}}$")

tt(as.data.frame(table_2),
   notes = "\\footnotesize \\textbf{Notes:} Table reports estimates from unweighted estimates. Weighted estimated can be found in Table ") |>
  group_tt(
    j = list("Binary" = 2:3,
             "Integer" = 4:5,
             "Integer" = 6:7)
  ) |>
  group_tt(
    j = list("Policy Support" = 2:5,
             "Trust" = 6:7)
    ) |>
  group_tt(
    j = list("$H_0:~\\beta_{\\text{Accountability}}-\\beta_{X}=0$" = 2:7)
  ) |>
  style_tt(j = 2:7, align = "c") |> 
  style_tt(j = 1, align = "r") |> 
  save_tt(here("output", "tables", "table-2.tex"),
          overwrite = TRUE)

# Figures 3 (and A6 and A7) -----------------------------------------------

survey2$policy_support_bin_int <- as.integer(survey2$policy_support_bin)
survey2$party_fct <- factor(survey2$party)
survey2$college_bin_fct <- factor(survey2$college_bin)
survey2$female_fct <- factor(survey2$female)
survey2$stock_fct <- factor(survey2$stock)


## Nationalism ------------------------------------------------------------

pdf(here('output', 'figures', 'figure-3a.pdf'),
    height = 4, width = 5, onefile=FALSE)
interflex(estimator = 'binning',
          Y = 'policy_support_bin_int',
          D = 'treat',
          base = 0,
          X = 'nat_total_int',
          Z = c('party_fct',
                'college_bin_fct', 'female_fct', 
                'age', 'stock_fct', 'income_int'),
          na.rm = T,
          vcov.type = 'robust',
          treat.type = 'discrete',
          nbins = 3,
          Xdistr = 'hist',
          theme.bw = TRUE,
          ylab = 'Marginal Effect of China Treatment',
          xlab = 'National Superiority Index',
          cex.lab = .85,
          cex.axis = .75,
          data = survey2[survey2$attentive == 1, ])
dev.off()

pdf(here('output', 'figures', 'figure-A6a.pdf'),
    height = 4, width = 5, onefile=FALSE)
interflex(estimator = 'binning',
          Y = 'policy_support_bin_int',
          D = 'treat',
          base = 0,
          X = 'nat_total_int',
          Z = c('party_fct',
                'college_bin_fct', 'female_fct', 
                'age', 'stock_fct', 'income_int'),
          weights = 'w_ebal_attentive',
          na.rm = T,
          vcov.type = 'robust',
          treat.type = 'discrete',
          nbins = 3,
          Xdistr = 'hist',
          theme.bw = TRUE,
          ylab = 'Marginal Effect of China Treatment',
          xlab = 'National Superiority Index',
          cex.lab = .85,
          cex.axis = .75,
          data = survey2[survey2$attentive == 1, ])
dev.off()


pdf(here('output', 'figures', 'figure-A7a.pdf'),
    height = 4, width = 5, onefile=FALSE)
interflex(estimator = 'kernel',
          Y = 'policy_support_bin_int',
          D = 'treat',
          base = 0,
          X = 'nat_total_int',
          Z = c('party_fct',
                'college_bin_fct', 'female_fct', 
                'age', 'stock_fct', 'income_int'),
          na.rm = T,
          vcov.type = 'robust',
          treat.type = 'discrete',
          Xdistr = 'hist',
          theme.bw = TRUE,
          ylab = 'Marginal Effect of China Treatment',
          xlab = 'National Superiority Index',
          cex.lab = .85,
          cex.axis = .75,
          data = survey2[survey2$attentive == 1, ])
dev.off()


## Active Engagement ------------------------------------------------------

pdf(here('output', 'figures', 'figure-3b.pdf'),
    height = 4, width = 5, onefile=FALSE)
interflex(estimator = 'binning',
          Y = 'policy_support_bin_int',
          D = 'treat',
          base = 0,
          X = 'internationalism_total_int',
          Z = c('party_fct',
                'college_bin_fct', 'female_fct', 
                'age', 'stock_fct', 'income_int'),
          na.rm = T,
          vcov.type = 'robust',
          treat.type = 'discrete',
          nbins = 3,
          Xdistr = 'hist',
          theme.bw = TRUE,
          ylab = 'Marginal Effect of China Treatment',
          xlab = 'Active Engagement Index',
          cex.lab = .85,
          cex.axis = .75,
          data = survey2[survey2$attentive == 1, ])
dev.off()

pdf(here('output', 'figures', 'figure-A6b.pdf'),
    height = 4, width = 5, onefile=FALSE)
interflex(estimator = 'binning',
          Y = 'policy_support_bin_int',
          D = 'treat',
          base = 0,
          X = 'internationalism_total_int',
          Z = c('party_fct',
                'college_bin_fct', 'female_fct', 
                'age', 'stock_fct', 'income_int'),
          weights = 'w_ebal_attentive',
          na.rm = T,
          vcov.type = 'robust',
          treat.type = 'discrete',
          nbins = 3,
          Xdistr = 'hist',
          theme.bw = TRUE,
          ylab = 'Marginal Effect of China Treatment',
          xlab = 'Active Engagement Index',
          cex.lab = .85,
          cex.axis = .75,
          data = survey2[survey2$attentive == 1, ])
dev.off()

pdf(here('output', 'figures', 'figure-A7b.pdf'),
    height = 4, width = 5, onefile=FALSE)
interflex(estimator = 'kernel',
          Y = 'policy_support_bin_int',
          D = 'treat',
          base = 0,
          X = 'internationalism_total_int',
          Z = c('party_fct',
                'college_bin_fct', 'female_fct', 
                'age', 'stock_fct', 'income_int'),
          na.rm = T,
          vcov.type = 'robust',
          treat.type = 'discrete',
          Xdistr = 'hist',
          theme.bw = TRUE,
          ylab = 'Marginal Effect of China Treatment',
          xlab = 'Active Engagement Index',
          cex.lab = .85,
          cex.axis = .75,
          data = survey2[survey2$attentive == 1, ])
dev.off()

# Figure 4 ----------------------------------------------------------------

m_zero_sum <- lm_robust(policy_support_bin ~ treat * zero_sum
                        + party 
                        + female
                        + age
                        + income_int
                        + stock 
                        + college_bin, 
                        se_type = 'HC2',
                        survey2[survey2$attentive == 1, ])

m_zero_sum_weighted <- lm_robust(policy_support_bin ~ treat * zero_sum
                        + party 
                        + female
                        + age
                        + income_int
                        + stock 
                        + college_bin, 
                        se_type = 'HC2',
                        weights = w_ebal_attentive,
                        survey2[survey2$attentive == 1, ])

plotData <- plot_slopes(m_zero_sum,
                        variables = 'treat',
                        condition = 'zero_sum',
                        draw = FALSE)

plotData_weighted <- plot_slopes(m_zero_sum_weighted,
                                 variables = 'treat',
                                 condition = 'zero_sum',
                                 draw = FALSE)

plotData$zs_label <- ifelse(plotData$zero_sum == 1, 'Yes', 'No')
plotData$Weighted <- "No"

plotData_weighted$zs_label <- ifelse(plotData_weighted$zero_sum == 1, 'Yes', 'No')
plotData_weighted$Weighted <- "Yes"

cols <- c("estimate", "conf.high", "conf.low", "std.error", "zs_label", "Weighted")

ggplot(rbind(plotData[, cols], plotData_weighted[, cols]), 
       aes(y = zs_label,
           x = estimate,
           color = Weighted,
           shape = Weighted)) +
  geom_vline(xintercept = 0,
             size = .75, 
             linetype = 'dashed',
             color = 'red') +
  geom_point(size = 4.5,
             position=position_dodge(width=-.5)) +
  geom_errorbar(aes(xmin = conf.high, xmax = conf.low),
                width = .15, 
                size = .5,
                position=position_dodge(width=-.5)) +
  geom_errorbar(aes(xmin = estimate + (std.error * 1.64),
                    xmax = estimate - (std.error * 1.64)),
                width = 0, size = 1.5,
                position=position_dodge(width=-.5)) +
  labs(x = 'Marginal Effect of Treatment',
       y = 'Zero Sum?') +
  scale_color_manual(values = c("black", "darkgrey")) +
  theme_bw(base_size = 15) +
  theme(legend.position = "right")

ggsave(here('output', 'figures', 'figure-4.pdf'),
       height = 3, width = 6.5, units = 'in')

# Table A1 -------------------------------------------

cols <- c('outcome_policy_support', 'outcome_policy_int', 'outcome_trust',
          'democrat', 'republican', 
          'college', 'female', 
          'stock', 'age', 'hhi', 'nat_total_int')

table <- survey1[attentive == 1, 
                 lapply(.SD, function (x) {paste0(format(round(mean(x), 3), nsmall = 3), 
                                                  ' (', 
                                                  format(round(sd(x), 3), nsmall = 3), 
                                                  ')' )
                 } ), 
                 .SDcols = cols, by = fed_treatment_policy]


table <- t(table)
colnames(table) <- c('Control', 'Risk', 'Accountability', 'Moral Hazard')

rownames(table) <- c('Variable', 
                     'Support (Binary)',
                     'Support (Integer)',
                     'Trust',
                     'Democrat',
                     'Republican',
                     'College',
                     'Female',
                     'Stock',
                     'Age',
                     'Income',
                     'Nationalism')

printTable <- xtable(table[-1, ], align = c('lcccc'),
                     caption = 'Policy Treatments, Lucid Sample',
                     label = 'table:desc-policy')

print(printTable, 
      include.rownames = TRUE,
      include.colnames = TRUE,
      booktabs = TRUE,
      file = here('output', 'tables', 'table-A1.tex'))

# Table A2 -----
table <- survey1[attentive ==1 , 
                 lapply(.SD, function (x) {paste0(format(round(mean(x), 3), nsmall = 3), 
                                                  ' (', 
                                                  format(round(sd(x), 3), nsmall = 3), 
                                                  ')' )
                 } ), 
                 .SDcols = cols, by = fed_treatment_china]


table <- t(table)
colnames(table) <- c('Control', 'China')

rownames(table) <- c('Variable', 
                     'Support (Binary)',
                     'Support (Integer)',
                     'Trust',
                     'Democrat',
                     'Republican',
                     'College',
                     'Female',
                     'Stock',
                     'Age',
                     'Income',
                     'Nationalism')

printTable <- xtable(table[-1, ], align = c('lcc'),
                     caption = 'China Treatments, Lucid Sample',
                     label = 'table:desc-china')

print(printTable, 
      include.rownames = TRUE,
      include.colnames = TRUE,
      booktabs = TRUE,
      file = here('output', 'tables', 'table-A2.tex'))

# Table A3 ------
cols <- c('policy_support_bin', 'policy_support_int', 'trust',
          'Republican', 'Democrat', 'Independent',
          'college_bin', 'female', 
          'age', 'stock', 'income_int', 
          'nat_total_int', 'internationalism_total_int', 'zero_sum')

table <- survey2[attentive == 1, lapply(.SD, function (x) {
  paste0(format(round(mean(x), 3), nsmall = 3), 
         ' (', 
         format(round(sd(x), 3), nsmall = 3), 
         ')' )
} ), .SDcols = cols, by = treat]


table <- t(table)
colnames(table) <- c('Control', 'China')

rownames(table) <- c('Variable', 
                     'Support (Binary)',
                     'Support (Integer)',
                     'Trust',
                     'Republican',
                     'Democrat',
                     'Independent',
                     'College',
                     'Female',
                     'Age',
                     'Stock',
                     'Income',
                     'Nationalism',
                     'Internationalism',
                     'Zero Sum')

printTable <- xtable(table[-1, ], align = c('lcc'),
                     caption = 'Descriptive Statistics, Prolific Sample',
                     label = 'table:desc-survey2')

print(printTable, 
      include.rownames = TRUE,
      include.colnames = TRUE,
      booktabs = TRUE,
      file = here('output', 'tables', 'table-A3.tex'))

rm(table, printTable, cols)

# Figure A1 -----

pdf(here('output', 'figures', 'figure-A1.pdf'),
    width = 5.5, height = 4)
hist(survey1$nat_total_int,
     xlab = 'Nationalism Index',
     main = '')
dev.off()

# Figure A2 ------

pdf(here('output', 'figures', 'figure-A2a.pdf'),
    width = 5.5, height = 4)
hist(survey2$nat_total_int,
     xlab = 'Nationalism Index',
     main = '')
dev.off()

pdf(here('output', 'figures', 'figure-A2b.pdf'),
    width = 5.5, height = 4)
hist(survey2$internationalism_total_int,
     xlab = 'Active Engagement Index',
     main = '')
dev.off()

# Tables A4 -----------------------------------------------

coefNames <- list(
  'fed_treatment_policyfed_treat_dem' = 'Accountability',
  'fed_treatment_policyfed_treat_mh' = 'Moral Hazard',
  'fed_treatment_policyfed_treat_risk' = 'Risk',
  'fed_treatment_china' = 'China',
  'democrat' = 'Democrat',
  'republican' = 'Republican',
  'college' = 'College',
  'female' = 'Female',
  'stock' = 'Own Stock',
  'age' = 'Age',
  'hhi' = 'Income',
  '(Intercept)' = 'Intercept'
)

texreg(models_svy1,
       custom.coef.map = coefNames,
       stars = c(0.01, 0.05, .1),
       digits = 2,
       include.ci = FALSE,
       custom.header = list("Policy Support" = 1:4, "Trust" = 5:6),
       custom.model.names = paste("(", 1:6, ")"),
       custom.gof.rows = list("Outcome" = list("(0,1)", "(0,1)",
                                               "(1-6)", "(1-6)",
                                               "(0-10)", "(0-10)"),
                              "Weighted" = list("No", "Yes",
                                                "No", "Yes",
                                                "No", "Yes")),
       custom.note = '\\textbf{Note:}~Robust standard errors in parentheses.',
       caption = 'Main Results (Lucid)',
       label = 'table:appendix_main_results_fed_svy1',
       booktabs = TRUE,
       dcolumn = TRUE,
       use.packages = FALSE,
       file = here('output', 'tables', 'table-A4.tex'))

coefNames <- list(
  'treat' = 'China',
  'partyIndependent' = 'Independent',
  'partyRepublican' = 'Republican',
  'college_binTRUE' = 'College',
  'femaleTRUE' = 'Female',
  'stockTRUE' = 'Own Stock',
  'age' = 'Age',
  'income_int' = 'Income',
  '(Intercept)' = 'Intercept'
)

# Table A5 -------
coefNames <- list(
  'treat' = 'China',
  'partyIndependent' = 'Independent',
  'partyRepublican' = 'Republican',
  'college_binTRUE' = 'College',
  'femaleTRUE' = 'Female',
  'stockTRUE' = 'Own Stock',
  'age' = 'Age',
  'income_int' = 'Income',
  '(Intercept)' = 'Intercept'
)

texreg(models_svy2,
       custom.coef.map = coefNames,
       stars = c(0.01, 0.05, .1),
       digits = 2,
       include.ci = FALSE,
       custom.header = list("Policy Support" = 1:4, "Trust" = 5:6),
       custom.model.names = paste("(", 1:6, ")"),
       custom.gof.rows = list("Outcome" = list("(0,1)", "(0,1)",
                                               "(1-6)", "(1-6)",
                                               "(0-10)", "(0-10)"),
                              "Weighted" = list("No", "Yes",
                                                "No", "Yes",
                                                "No", "Yes")),
       custom.note = '\\textbf{Note:}~Robust standard errors in parentheses.',
       caption = 'Main Results (Prolific)',
       label = 'table:appendix_main_results_fed_svy2',
       booktabs = TRUE,
       dcolumn = TRUE,
       use.packages = FALSE,
       file = here('output', 'tables', 'table-A5.tex'))

# Figure A3 ---------------------------------------------------------

## Analysis ---------------------------------------------------------------

### Policy Support, Binary, Unweighted -----
svy1_policy_bin_atten <- lh_robust(
  outcome_policy_support ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi
  + attentive,
  data = survey1,
  linear_hypothesis = c('fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh',
                        'fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk'),
)

models_svy1 <- append(models_svy1, list("Support, Binary, Unweighted, w/ inattantive" = svy1_policy_bin_atten))

### Policy Support, Binary, Weighted -----
svy1_policy_bin_w_atten <- lh_robust(
  outcome_policy_support ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi
  + attentive,
  weights = w_ebal_full,
  data = survey1,
  linear_hypothesis = c('fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh',
                        'fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk'),
)

models_svy1 <- append(models_svy1, list("Support, Binary, Weighted, w/ inattantive" = svy1_policy_bin_w_atten))

### Policy Support, Integer, Unweighted -----
svy1_policy_int_atten <- lh_robust(
  outcome_policy_int ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi
  + attentive,
  data = survey1,
  linear_hypothesis = c('fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh',
                        'fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk'),
)

models_svy1 <- append(models_svy1, list("Support, Integer, Unweighted, w/ inattantive" = svy1_policy_int_atten))

### Policy Support, Integer, Weighted -----
svy1_policy_int_w_atten <- lh_robust(
  outcome_policy_int ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi
  + attentive,
  data = survey1,
  weights = w_ebal_full,
  linear_hypothesis = c('fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh',
                        'fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk'),
)

models_svy1 <- append(models_svy1, list("Support, Integer, Weighted, w/ inattantive" = svy1_policy_int_w_atten))

### Trust, Integer, Unweighted -----
svy1_trust_atten <- lh_robust(
  outcome_trust ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi
  + attentive,
  data = survey1,
  linear_hypothesis = c('fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh',
                        'fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk'),
)

models_svy1 <- append(models_svy1, list("Trust, Integer, Unweighted, w/ inattantive" = svy1_trust_atten))

### Trust, Integer, Weighted -----
svy1_trust_w_atten <- lh_robust(
  outcome_trust ~ fed_treatment_policy
  + fed_treatment_china 
  + democrat 
  + republican 
  + college 
  + female 
  + stock 
  + age 
  + hhi
  + attentive,
  weights = w_ebal_full,
  data = survey1,
  linear_hypothesis = c('fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh',
                        'fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk'),
)

models_svy1 <- append(models_svy1, list("Trust, Integer, Weighted, w/ inattantive" = svy1_trust_w_atten))


## Plot Figure A3 --------------------------------------------------------

params_svy1 <- c('Accountability' = 'fed_treatment_policyfed_treat_dem',
                 'Moral Hazard' = 'fed_treatment_policyfed_treat_mh',
                 'Risk' = 'fed_treatment_policyfed_treat_risk',
                 'China' = 'fed_treatment_china')

plotData <- data.table()
results <- data.table()

for (m in 1:length(models_svy1)) {
  results <- rbind(results,
                   cbind("model" = names(models_svy1)[m],  
                         get_cis(model = models_svy1[[m]]$lm_robust,
                                 params = params_svy1,
                                 levels = c(.95, .90))
                   ) # end cbind
  ) # end rbind
}

results$outcome <- sapply(strsplit(results$model, ", "),
                          "[[", 1)
results$outcome_measure <- sapply(strsplit(results$model, ", "),
                                  "[[", 2)
results$weight <- sapply(strsplit(results$model, ", "),
                         "[[", 3)

results$attentive <- ifelse(grepl("inattantive", results$model,fixed = TRUE),
                            "Yes", "No")

results$var_name <- factor(results$var_name, levels = c('China',
                                                        'Risk',
                                                        'Moral Hazard',
                                                        'Accountability'))

results$outcome_label <- paste0(results$outcome, " (", results$outcome_measure, ")")

ggplot(results[results$attentive == "Yes", ], aes(y = var_name, x = est, 
                    color = weight,
                    shape = weight)) +
  geom_vline(xintercept = 0, color = 'red', linetype = 'dashed', size = .5) +
  geom_point(size = 3.5,
             position=position_dodge(width=-1)) +
  geom_errorbar(aes(xmin = `2.5 %`, xmax = `97.5 %`),
                #width = c(rep(.2, 6), rep(.4, 4)), 
                width = 0,
                size = .66,
                position=position_dodge(width=-1)
  ) +
  geom_errorbar(aes(xmin = `5 %`, xmax = `95 %`),
                width = 0, 
                size = 1.45,
                position=position_dodge(width=-1)
  ) +
  scale_color_manual(name = "Weighted?",
                     values = c('black', 'gray')) +
  scale_shape_manual(name = "Weighted?",
                     values = c('circle', 'triangle')) +
  theme_bw(base_size = 12) +
  labs(y = element_blank(),
       x = element_blank(),
       color = 'Sample') +
  theme(legend.position = c(0.75, 0.2)) +
  facet_wrap(~outcome_label, 
             ncol = 2,
             scales = "free_x")

ggsave(here("output", "figures", "figure-A3.png"),
       width = 6.5, height = 5, dpi = 450, units = "in")

# Figure A4 -------
covars <- c('democrat', 'republican', 
                  'college', 'female', 
                  'stock', 'age', 'hhi' , 'fed_treatment_china')

newData <- survey1[ , lapply(.SD, mean), .SDcols = covars]

preds <- data.frame()

for (t in unique(survey1$fed_treatment_policy)) {
  newData$fed_treatment_policy <- t
  temp <- 
    preds <- rbind(preds, 
                   cbind(t, 
                         predict(models_svy1[[1]]$lm_robust, newData, 
                                 interval = "confidence")$fit
                   ) #cbind end
    ) #rbind end
}

newData$fed_treatment_china <- 1
newData$fed_treatment_policy <- 'fed_treat_lend'

preds <- rbind(preds, cbind(t='fed_treat_china', 
                            predict(models_svy1[[1]]$lm_robust, newData, 
                                    interval = 'confidence')$fit))

preds$treats <- c('Moral\nHazard', 'Control', 'Risk',
                  'Accountability', 
                  'China')

preds$treats <- factor(preds$treats, levels = c('Control', 'Risk',
                                                'Moral\nHazard', 
                                                'Accountability', 
                                                'China'))

preds[,2:4] <- apply(preds[,2:4], 2, as.numeric)

ggplot(preds, aes(x = treats, y = fit)) +
  geom_hline(yintercept = .5, color = 'red', linetype = 'dashed') +
  geom_bar(stat = 'identity', width = .5,
           fill = 'gray', color = 'black') +
  geom_errorbar(aes(ymin = lwr, ymax = upr),
                width = .2) +
  scale_y_continuous(limits = c(0, .593),
                     labels = scales::percent_format(1L),
                     breaks = seq(0, .6, .1),
                     expand = expansion(mult = c(0.0015, 0.04))) +
  theme_bw(base_size = 14) +
  labs(y = '% Support',
       x = element_blank())
rm(t)

ggsave(here('output', 'figures', 'figure-A4.pdf'),
       height = 4, width = 6, units = 'in')

# Table A6 -----------

## Lending ----------------------------------------------------------------

model_policy_bin_lin <- lm_lin(
  outcome_policy_support ~ fed_treatment_policy,
  covariates = ~ democrat+republican+college+female+stock+age+hhi,
  se_type = "HC2",
  data = survey1[survey1$attentive == 1,]
)

model_policy_int_lin <- lm_lin(
  outcome_policy_int ~ fed_treatment_policy,
  covariates = ~ democrat+republican+college+female+stock+age+hhi,
  se_type = "HC2",
  #weights = w_ebal_attentive,
  data = survey1[survey1$attentive == 1,]
)

model_trust_int_lin <- lm_lin(
  outcome_trust ~ fed_treatment_policy,
  covariates = ~ democrat+republican+college+female+stock+age+hhi,
  se_type = "HC2",
  #weights = w_ebal_attentive,
  data = survey1[survey1$attentive == 1,]
)

## China ------------------------------------------------------------------

model_policy_bin_china_lin <- lm_lin(
  outcome_policy_support ~ fed_treatment_china,
  covariates = ~ democrat+republican+college+female+stock+age+hhi,
  se_type = "HC2",
  data = as.data.frame(survey1[survey1$attentive == 1,])
)

model_policy_int_china_lin <- lm_lin(
  outcome_policy_int ~ fed_treatment_china,
  covariates = ~ democrat+republican+college+female+stock+age+hhi,
  se_type = "HC2",
  data = as.data.frame(survey1[survey1$attentive == 1,])
)

model_trust_int_china_lin <- lm_lin(
  outcome_trust ~ fed_treatment_china,
  covariates = ~ democrat+republican+college+female+stock+age+hhi,
  se_type = "HC2",
  data = as.data.frame(survey1[survey1$attentive == 1,])
)

## Make table -------------------------------------------------------------

coefNames <- list(
  'fed_treatment_policyfed_treat_dem' = 'Accountability',
  'fed_treatment_policyfed_treat_mh' = 'Moral Hazard',
  'fed_treatment_policyfed_treat_risk' = 'Risk',
  'fed_treatment_china' = 'China'
)

texreg(list(model_policy_bin_lin,
            model_policy_int_lin,
            model_trust_int_lin,
            model_policy_bin_china_lin,
            model_policy_int_china_lin,
            model_trust_int_china_lin),
       custom.coef.map = coefNames,
       stars = c(0.01, 0.05, .1),
       digits = 3,
       include.ci = FALSE,
       custom.header = list("Policy" = 1:3, "China" = 4:6),
       custom.model.names = rep(c('Support', 'Support', 'Trust'), 2),
       custom.note = '\\textbf{Note:} Robust standard errors in parentheses.',
       caption = 'Results from estimator proposed by \\cite{linAgnosticNotesRegression2013}',
       label = 'table:appendix_main_lin_2013',
       booktabs = TRUE,
       dcolumn = TRUE,
       use.packages = FALSE,
       file = here('output', 'tables', 'table-A6.tex'))

# Figure A5 -----

results_a5 <- data.table()

# Finish
params_svy1_a5 <- c('Moral Hazard' = "fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_mh",
                    'Risk' = "fed_treatment_policyfed_treat_dem = fed_treatment_policyfed_treat_risk")

for (m in 1:length(models_svy1)) {
  results_a5 <- rbind(results_a5,
                   cbind("model" = names(models_svy1)[m],  
                         get_cis(model = models_svy1[[m]]$lh,
                                 params = params_svy1_a5,
                                 levels = c(.95, .90),
                                 lh = TRUE)
                   ) # end cbind
  ) # end rbind
}

results_a5$outcome <- sapply(strsplit(results_a5$model, ", "),
                          "[[", 1)
results_a5$outcome_measure <- sapply(strsplit(results_a5$model, ", "),
                                  "[[", 2)
results_a5$weight <- sapply(strsplit(results_a5$model, ", "),
                         "[[", 3)

results_a5$attentive <- ifelse(grepl("inattantive", results_a5$model,fixed = TRUE),
                            "Yes", "No")

results_a5$var_name <- factor(results_a5$var_name, levels = c('Risk',
                                                        'Moral Hazard'))

results_a5$outcome_label <- paste0(results_a5$outcome, " (", results_a5$outcome_measure, ")")

# Plot differences accross all models

# Trust
p1 <- ggplot(results_a5[results_a5$outcome == "Trust", ],
       aes(x = est,
           y = model,
           shape = attentive,
           color = weight)) +
  geom_vline(xintercept = 0, color = 'red', linetype = 'dashed', size = .5) +
  geom_point(size = 3.5) +
  geom_errorbar(aes(xmin = (est - (se * 1.96)), 
                    xmax = (est + (se * 1.96))),
                width = 0,
                size = .66,
                position=position_dodge(width=-1)
  ) +
  geom_errorbar(aes(xmin = (est - (se * 1.64)), 
                    xmax = (est + (se * 1.64))),
                width = 0, 
                size = 1.45,
                position=position_dodge(width=-1)
  ) +
  labs(title = "Trust",
       shape = "Inc. inattentive respondents?",
       color = "Weighted?",
       x = element_blank(),
       y = element_blank()) +
  theme_bw(base_size = 11) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        legend.title.position = "top") +
  facet_wrap(~coef)

# Support, Integer
p2 <- ggplot(results_a5[results_a5$outcome == "Support" & results_a5$outcome_measure == "Integer", ],
       aes(x = est,
           y = model,
           shape = attentive,
           color = weight)) +
  geom_vline(xintercept = 0, color = 'red', linetype = 'dashed', size = .5) +
  geom_point(size = 3.5) +
  geom_errorbar(aes(xmin = (est - (se * 1.96)), 
                    xmax = (est + (se * 1.96))),
                width = 0,
                size = .66,
                position=position_dodge(width=-1)
  ) +
  geom_errorbar(aes(xmin = (est - (se * 1.64)), 
                    xmax = (est + (se * 1.64))),
                width = 0, 
                size = 1.45,
                position=position_dodge(width=-1)
  ) +
  labs(title = "Support, Integer",
       shape = "Inc. inattentive respondents?",
       color = "Weighted?",
       x = element_blank(),
       y = element_blank()) +
  theme_bw(base_size = 11) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        legend.title.position = "top") +
  facet_wrap(~coef)

# Support, Binary
p3 <- ggplot(results_a5[results_a5$outcome == "Support" & results_a5$outcome_measure == "Binary"],
       aes(x = est,
           y = model,
           shape = attentive,
           color = weight)) +
  geom_vline(xintercept = 0, color = 'red', linetype = 'dashed', size = .5) +
  geom_point(size = 3.5) +
  geom_errorbar(aes(xmin = (est - (se * 1.96)), 
                    xmax = (est + (se * 1.96))),
                width = 0,
                size = .66,
                position=position_dodge(width=-1)
  ) +
  geom_errorbar(aes(xmin = (est - (se * 1.64)), 
                    xmax = (est + (se * 1.64))),
                width = 0, 
                size = 1.45,
                position=position_dodge(width=-1)
  ) +
  labs(title = "Support, Binary",
       shape = "Inc. inattentive respondents?",
       color = "Weighted?",
       x = element_blank(),
       y = element_blank()) +
  theme_bw(base_size = 11) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        legend.title.position = "top") +
  facet_wrap(~coef)

ggpubr::ggarrange(p1, p2, p3, ncol = 1,
                  common.legend = TRUE,
                  legend = "bottom")

ggsave(here("output", "figures", "figure-A5.png"),
       width = 7, height = 9, dpi = 350)


